Smoothly-varying hopping rates in driven flow with exclusion 
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We consider the one-dimensional totally asymmetric simple exclusion process (TASEP) with 
position-dependent hopping rates. The problem is solved, in a mean field/adiabatic approxima- 
tion, for a general (smooth) form of spatial rate variation. Numerical simulations of systems with 
hopping rates varying linearly against position (constant rate gradient), for both periodic and open 
boundary conditions, provide detailed confirmation of theoretical predictions, concerning steady- 
state average density profiles and currents, as well as open-system phase boundaries, to excellent 
numerical accuracy. 

PACS numbers: 05.40.-a,02.50.-r,05.70.Fh 



I. INTRODUCTION 

In this paper we investigate the one-dimensional totally 
asymmetric simple exclusion process (TASEP) fll|, in the 
presence of non-uniform hopping rates. The TASEP is 
a biased diffusion process forparticles with hard-core re- 
pulsion (excluded volume) Notwithstanding the 
simplicity of formulation of its basic rules, this model can 
exhibit a wealth of non-trivial properties, and is consid- 
ered a paradigm in the field of non-equilibrium phenom- 
ena. Quenched random inhomogeneities in the TASEP 
have been extensively considered earlier [5l-[l2j. In con- 
trast, the case of deterministically- varying, position- 
depende nt p hysical parameters has received less atten- 
tion 

The TASEP and its generalizations have been applied 
to a broad range of non-equilibrium physical contexts, 
from macroscopic ones such as highway traffic [l5| to mi- 
croscopic ones, including sequence alignment in computa- 
tional biology [l6l | and current shot noise in quantum-dot 
chains [T7| . Situations may arise where monotonic spa- 
tial variations in an associated parameter can be relevant 
(such as gradients in the first case, and the "gap-cost", 
or an applied electric field, for the latter two cases). By 
contrast, the effects of, e.g., temperature gradients on 
the equilibrium [18] and transport properties of spin 
systems have been studied in detail; the same applies to 
concentration gradients in percolation (20l - l22| . One typ- 
ically gets a picture of spatial phase separation, in which 
a high-temperature (or low-concentration) disordered re- 
gion connects to a low-temperature (high-concentration) 
ordered one via an interface, whose features (e.g., width) 
scale in a non-trivial way with the inhomogeneity pa- 
rameters. More recently, experimental progress in cold- 
atom trapping [23] has been one motivation behind the 
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theoretical study of (pseudo)-spin systems in trapping 
potentials such as roagnetic fields with a wedge-like or 
parabolic profile |24l - l26l ]. 

We consider the problem of flow with exclusion, for 
which the time evolution of the 1-1-1 dimensional TASEP 
is the fundamental discrete model. The particle number 
n£ at lattice site i can be or 1, and the forward hop- 
ping of particles is only to an empty adjacent site. The 
current across the bond from £ to £ -|- 1 depends also on 
the stochastic attempt rate, pi, associated to it and is 
thus given by Jij+i — peng{l — n^+i) . For the usual ho- 
mogeneous case of p£ = p, in numerical simulations one 
can effectively make p — I, provided that the inherent 
stochasticity of the process is kept, via e.g. random se- 
lection of site occupation update [llj. This amounts to 
a trivial renormalization of the time scale. 

Here, we consider a position-dependent hopping rate 
(which cannot thus be simply renormalized away) . By us- 
ing periodic or open boundary conditions, with assorted 
overall densities in the former case, and injection/ejection 
rates in the latter, we investigate the consequent effects 
upon the associated particle density profiles and currents. 

To begin with we give the generic dynamic mean field 
theory for arbitrary "slow" space-dependence of the hop- 
ping rate. We then turn, for more specific results, to the 
steady state in the case of a linear dependence of pi on 
position (uniform gradient). It is remarkable that, from 
the combination of the mean field approach with an adia- 
batic approximation (to be described below) , many accu- 
rate results are obtained, including some such as current, 
and open-system phase boundaries, which appear to be 
exact in the large-system limit. 

Section |lT] below gives the mean-field/ adiabatic the- 
ory. In Section Hill we investigate the TASEP with peri- 
odic boundary conditions; in Sec. llVi we examine open- 
boundary TASEP systems in the following phases: (a) 
maximal-current, (b) low-density, (c) high-density, and 
(d) on the coexistence line. Finally, in Sec. |Vl conclud- 
ing remarks are made. 
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II. PRELIMINARIES AND BASIC THEORY 
A. Preliminaries 

We start by imposing periodic boundary conditions 
(PBC) for the TASEP at the ends of the chain, thus the 
total number of particles is fixed. For a uniform system 
in the steady state, the local average density at all sites 
coincides with the position-averaged particle density (p) 
(also to be denoted below by p, wherever no chance of a 
misunderstanding arises) . 

Although this is a discrete model, we denote positions 
along the lattice by a continuous variable x, thus (with 
the lattice parameter being of unit length), the bond la- 
belled by X connects sites x — ^ and x + ^ . The use of 
a continuum description is consistent with our emphasis 
throughout the paper on results applying in the infinite- 
system limit. 

We consider a linearly-varying hopping rate; although 
the theory developed in Subsection IIIBI below applies 
to a general position dependence (provided some rather 
general smoothness assumptions are valid), this constant- 
gradient case will be our choice of concrete application in 
the subsequent sections. For a system of size L, we take 



p{x) =po + 9- 



L L 
2 - - 2 



(1) 



where 9 denotes the intensity of the hopping-rate gradi- 
ent; we keep pq = 1/2 henceforth. 

The effect of the hopping-rate gradient, given by 
Eq. ([!]), on local densities is rather remarkable, as il- 
lustrated in Figure [1] . 

A schematic interpretation of the profile shape dis- 
played in Figure [T] can be provided as follows, using 
ideas from previous treatments of the quenched random- 
bond version of the TASEP d M^E^- For the TASEP 
with uniform rates p, it is known [lj-|j| that, for currents 
greater or less than Jc (p) = j the steady state phases are 
characterized by density profiles which are either: mono- 
tonically decreasing, {p{x)) — ^ = —q tanq(x — xo) (high- 
current phase) or monotonically increasing, {p{x)) — ^ = 
k tanhfc(a; — Xq) (kink-like, low-current phase). Here, q 
and k are characteristic inverse lengths such that q^ = 
-fc2 = (J _ JJ/4p [1, [13, [13, where J is the steady- 
state current; the profile forms result from the fact that 
J is constant throughout the system. This latter fact 
has strong bearing on the local shape of density pro- 
files in the quenched random-bond case: in regions with 
weak (strong) bonds, i.e. bonds with low (high) hop- 
ping probability (Ps), J can be larger (smaller) than 
the local critical current Jc{pw) (JciPs)), in which case 
the profile is of high-current (low-current) type. With 
9 > in Eq. ([T]), the features shown in Figure [1] ap- 
pear roughly consistent with the theoretical framework 
just sketched. However, we shall see from the full treat- 
ment developed in Subsection IH Bl below that, although 
the concepts of high- and low-current phases still per- 
sist here, their effects are strongly modified by factors 
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Figure 1. (Color online) Points give steady-state density 
profile for TASEP with hopping-rate gradient, for periodic 
boundary conditions, (p) = 1/2, lattice size L = 256, and 
9 — 0.2 (see Eq. Full lines are fits to forms inspired by 

the theory of randomly-disordered systems (see text). 



specific to the present case. In particular, the separa- 
tion in space of the two phases is actually very close to 
the left boundary in Figure [TJ not where the tan and 
tanh functions join in the fit shown in that same Figure. 
This is because the actual profiles involve tan and tanh 
functions with spatially varying "envelope" factors (see 
Subsection IH B[) . Many new features will be seen to arise 
from the "registration" in space of the envelope, i.e. its 
position in the system; as we shall see, the location of the 
envelope relative to the region of weakest bonds is set by 
the current. 

From the conjunction of PBC with the form of p{x) 
given in Eq. ([1]), one sees that particles find a hopping- 
rate discontinuity of amplitude —9 as they jump across 
the chain's endpoint. Although, from elementary consid- 
erations, PBC impose continuity of p across the gap, it is 
important to emphasize that the kink-like profile seen in 
Figure [1] is not an artifact brought about by the discon- 
tinuity just mentioned. As we shall see in the following, 
kinks may (or may not) be present with PBC. Their exis- 
tence, or lack thereof, depends on combinations of p and 
9 according to mechanisms described by our theory. 

One should note that, if the sign of 9 is reversed in 
Eq. the plot of (p) — 1/2 versus x simply gets point- 
refiected relative to the origin. 

The steady-state currents in the type of system stud- 
ied here also differ markedly from their uniform counter- 
parts. We recall that, for the latter with PBC, the rela- 
tionship between current Jq = J {9 = 0^p,p) and density 
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Figure 2. (Color online) Steady-state currents J against gra- 
dient intensity 9 for a system with L — 2048 and PBC, for 
densities as shown. Each point is an average taken over 100 
independent samples, each in its turn containing 1.2 x 10^ 
successive steady-state configurations. Error bars are smaller 
than symbol sizes. 



Jo =_pp(l - p) 



(2) 



where p is the uniform hopping rate. Eq. ^ is one ex- 
ample of relationships and quantities which mean field 
factorization gives exactly d, H^l , and whose generaliza- 
tion for non-uniform rates is also exactly given by the 
generalized mean field theory developed here, as we shall 
see. 

For now, we restrict ourselves to p < ^ . The question 
of whether or not the J — p diagram here displays the 
same symmetry, relative to p = i, as that for the uni- 
form case will be discussed later, with help of the theory 
developed in Subsection IIIBI . The effects on the cur- 
rent of a position-dependent p{x) given by Eq. JT]) are 
shown in Figure [51 for various densities, all of them not 
far removed from p — 1/2. It is seen that as 9 increases, 
the J — 9 relationship becomes independent of p for an 
increasingly broad range of densities, following the same 
nearly linear form as that of a system with p — 1/2. 
In other words, for fixed 9 a plateau develops around 
p = 1/2 in the J ~ p diagram, whose width increases 
with 9 . Again, a similar effect is seen in TASEP with 
quenched randomness @, IToj . 



B. Mean- field theory 

For uniform p, the Burgers equation 128143 1| , linearized 
via the Cole-Hopf transformation (3^ . l33l | gives the gen- 
eral time-dependent mean field solution, analogous to a 
superposition of moving solitons, corresponding to waves 
in the linearized system, with possibly complex wave- 
vectors. Real and imaginary wave vectors distinguish 
the two general soliton-like steady states which, because 
of particle conservation, are uniform-current ones. These 
steady states correspond to phases of maximal current 
(J ~ Jc = p/4) or low current (0 < J < Jc, and low 
or high density); the square of the wave vector is propor- 
tional to (J/ Jc — 1). For the special case of PBC, the two 
steady states become states of uniform density, while for 
open boundary conditions the steady state profiles are of 
tan and tanh form. 

For space-dependent p{x), the solution given below (for 
general time-dependence and then steady state) uses an 
adiabatic generalisation of constant-p ideas. 

We start from the continuity equation: 



dp 9 

with (using a mean field factorization) 

J{x) = p{x) p{x - 1-) (l-p{x + ^) 



(3) 



(4) 



defining a via p = ^(1 + 0-), Eq. ^ becomes, upon taking 
the continuum limit on Eq. 

'm=d-x^P^'^^^^'^ + d-x-'^^- 

Using the Cole-Hopf transformation (32l . H^j, we intro- 
duce the auxiliary variable u via a = dlnu/dx, in terms 
of which, after a partial integration with respect to x, 
Eq. turns into the linear form: 



(6) 



2--/(i).^p(.)|— 

where f{t) is the integration "constant". Writing u — 
X{x)T{t), one has 

whence 

r = exp(-ic.i + F(t)) with ^ = l/(t). (8) 



Putting 



{p^{x)f = 1 - 



p{x) 



and making the ansatz X = e'^'--^\ one gets 
c?7 



dx 



±Pu{x) , 



(9) 



(10) 
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provided d^i^(x)/dx <C (/itj(x))^ {adiabatic approxima- 
tion). In this limit X = e^'''-''^^\ with "fuj{x) = 
fj,ui{x) dx. Thus, 

XT = exp(±7^(x) - ^ujt + F{t)) . (11) 
The general solution for u(x,t) is 

u = ^A^cosh(7„(x)-7^(a^)) e-5-*+^^(*) , (12) 

where the A,^ and are arbitrary constants. Finally, in 
terms of (t(x, t): 

( .^ E,^^<^M'^(2;)sinh(7<^(a;) -7<^(a^)) e-^"* 

cr(x,rj = ;^ , 

cosh (7^(0;) - 7„(a„)) e-s-^* 

(13) 

in the mean field/adiabatic approximation. 
The following comments are in order: 

(i) If we take a single component in Eq. the e~2"* 
factor cancels and we are left with a steady state solution: 

cr(x) = ^i^(x)tanh(7„(a;) - 7^(a)) . (14) 

When the validity criterion for the adiabatic approxima- 
tion applies, this state is associated with the current 

J=ip(x)(l-M^(a;)2) = iw, (15) 

[using Eqs. (U)), ([H]), and (|14p ], which is constant as nec- 
essary for the steady state. 

(ii) In the t-dependent general form Eq. (jl3p . each sum 
evolves for long times into a single component, which is 
the one having the least w, corresponding to the steady 
state, i.e., = 4 J , by (i). 

(iii) At long but not infinite times the sums in Eq. (|13p 
are dominated by the terms with the smallest cj's. Then 
the denominator, whose logarithmic derivative gives tr, 
becomes a combination of the steady state component 
and a wave packet whose group velocity v{x) can be ob- 
tained by a straightforward adiabatic generalisation of 
standard procedures, using the analogue d^,^{x)/dx of 
the wave vector. The result is, generally, 

v[x)=±p{x)[l--^^' , (16) 

becoming v{x) = ±p{x)^{p{x) — 4 J) 5 for the kink dy- 
namics in the late-time approach to the steady state. 

In what follows we shall be mostly concerned with the 
steady state, so the following distinctions and details may 
be helpful. In Eq. (HI]), 

J,{x) = -^p{x) (17) 

acts like a local critical current, since the sign of J — Jc 
determines whether pi{x) there is real or imaginary and. 



consequently, whether the profile in Eq. involves a 
tanh or tan function. This is a generalization of the case 
with space-independent rate p, where Jc — p/4 is the 
maximal current, associated to flat or tan profiles, while 
low currents J < Jc exhibit tanh profiles. 

For the space-dependent p{x) the most important new 
features are the a;-dependence of Jc{x), the J-dependent 
location (xq) of the division between phases, and the oc- 
currence of the space-dependent amplitude function fj.{x) 
in the profile, Eq. (|14l) . Where it is necessary, to avoid 
confusion, we distinguish the possibilities by using, in 
place of p{x), the specific real functions k{x), q{x) de- 
fined by 

k{x) = fi{x) = (1 - U/p{x))i = (1 - J/Jc{x))-- , 

J < Jcix) ; 

q{x) ^l^i{x) = {UIp(x) - 1)^ = (J/Jc(.t) - 1)5 , 

J > Jc{x) . (18) 

Then 

a{x) = k{x) tanh(ii'(x) - K{a)), J < Jc{x) ; 
(t{x) = -q{x) tan{Q{x) - Q{b)), J > Jc{x) , (19) 

where K{x) = k{x) dx and Q{x) = q{x) dx. For 
x-dependent rates, the tan form can only apply in at 
most a very limited region (of size set by the weakest 
rates). This is because the tan function in a{x) diverges, 
violating the physical requirement on the local density, 
< Ij unless its argument Q{x) is limited to a range 
less than tt. So the tanh form will actually account for 
most of the profile. If the integration constant a is inside 
the system the change of sign of the argument of the tanh 
function at x = a corresponds to a kink there. For the 
tanh, k{x) acts like an envelope, and its crucial effects 
in distinguishing scenarios and phases partly relate to its 
registration, for which the tan part of the profile can play 
a dominant role. 



III. STEADY STATE WITH PBC 

A. Introduction 

For the non-uniform system with x € [— L/2,L/2], 
PBC impose the constraint on a{x): 

(t{-L/2) ^a{L/2) . (20) 

In addition to this, in order to fix arbitrary constants and 
determine the steady state current J and profile <j{x), 
we need also to specify the average density (p), in the 
equation 

2{{p)^\) = {a) = \ f'\{x)dx. (21) 
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1 
L 



With the mean field/adiabatic approximation this be- 
comes 

rL/2 

/ M(a^) tanh (7(2;) — 7(a)) dx = 

'-L/2 

= [\ncosMl{x)-l{a))t%^ . (22) 

Here we used /i = drjjdx^ and have reverted to non- 
specific notation, not distinguishing tanh or tan (nor cosh 
or cos) . We will later have to verify that the criterion for 
use of the adiabatic approximation is satisfied. 

From the general formulation above, /i and hence 7 
are related to the current J; it and the other parameter 
a (the kink position in the case of real 7(a)) are deter- 
mined in terms of (p) by Eqs. ((2T|) and ([22)) (for large 
systems, the kink position will be sharp when the adia- 
batic approximation is satisfied) . 



B. Rate gradient 

From now on we deal with the specific case of linearly- 
varying -pix) given in Eq. (P) . With PBC and > 0, 
one gets in the adiabatic approximation, with the help of 
Eqs. dHD, Cni), and (Uni): 



\i{x) 



4J 

p(x) _ 



X 



X 



(23) 



and 



^i{x) dx = [X{X + c)] 2 - ctanh 



X 



X + c\ 
K{X) , (24) 



where 



X — X ^ Xf) , 
a;o = (8J-l)^ 



4JL 



2^0 



L 
20 



-4 

L 

" 2 



- A 



(25) 



xq corresponds to the place where /i(x) vanishes, hence to 
the position of the apex of the envelope function ±|^(a;)|, 
i.e., where ^i{x) [and 7(2;) [ cross over between real and 
imaginary values k(x) or —i q(x) [ and K[x) or —i Q{x) ]. 
Subsequently explicit forms will be needed, particulary 
for 7 for the real case, and it will then be convenient to 
use both K and (real) K, where 

Kix) = K{X) , (26) 

with X ^ X{x) ^ X — xq, and where K is as in Eq. (IMl) . 
xq also corresponds to the place where J is equal to the 
local critical current; this plays a central role in the dis- 
cussion. For graphical illustrations, refer to Figure HI in 
subsection llll DI below. c is a characteristic length related 
to the rate gradient. A, the ratio of xq to —L/2, conve- 
niently distinguishes scenarios, and parametrizes analytic 
expressions, particulary in the L ^ 00 limit. 



C. Scenarios for steady state behavior 

We next discuss the character and location of steady 
state phases, and relationships to positions of the "enve- 
lope" and kinks. The generalised maximal current and 
low current phases of the system turn out to be described 
by two scenarios, I and II, as follows. 

For the rate gradient case with 6* > (0 < has dual 
character), the smallest p{x) is at the left-hand side edge, 
giving a severe bottleneck there. As we shall see in the 
following, this has the consequence that the current J 
adjusts itself in such a way that the apex position xq — 
—XL/2 turns out to be either: (I) near the left boundary, 
but still inside the system, or (II) to the left of the left 
boundary. These give, respectively: 

Scenario I: A < 1. Here the tan function applies near 
the left edge and its spatial extent Ax is limited by the 
condition AQ{x) ~ q{x) Ax < tt. Since q{x) is related 
to the difference J — Jc{x) of the steady state current J 
from its local critical value, this condition also limits J 
as well as the position, x = xq, where J — Jc{x) vanishes. 
Scenario II: A > 1. Here only the tanh function applies 
inside the system. 

The two scenarios become very evident in the "family" 
of profiles corresponding to all possible average densities 
(p), for PBC and a given 6 (see the numerical results in 
Figure [4]) . 

Scenario I corresponds to a common envelope (nearly 
parabolic in shape, see Eq. (j23p ) and applies for an in- 
termediate range of (p)'s (not very far from 1/2). It 
is consistent with a fixed position xq of the apex of 
the envelope, close to the left hand boundary. It is 
(through Eq. ([^5]) ) consistent with an observed constant 
((p) -independent) plateau current J, about {1/8){1 — 9). 
Near the left boundary there is a small region of tan pro- 
file, and everywhere else the profile approaches the tanh 
form (including the kink). 

Scenario II, applying for larger \{p) — 1/2|, has profiles 
not near a common envelope, corresponding to varying 
apex position; indeed, in this case xq is outside of the 
system (to the left of the left boundary) and the profile 
is entirely of tanh type. In this scenario the currents 
depend on (p). 

These scenarios and related phenomena can be quan- 
titatively explained using the mean field adiabatic for- 
mulation, except near the envelope apex if that lies in- 
side the system. This is because the apex is where fi{x) 
vanishes, i.e., where the adiabatic approximation fails 
utterly (see the validity criterion, below Eq. UO])). To 
the right of the apex, where J < Jc{x), the adiabatic 
approximation is valid for X > c^^^, so the adiabatic 
form an = k{x) tanh(_ftr(a;) — K(a)) applies; similarly, 
in the region to the left of the apex, J > Jc{x), and 
the adiabatic form = —q{x) tan((5(a;) — Q[h)) is valid 
for X < —c^/^ . Between these a (nonadiabatic) form 
ac cc {x + const. )^^ is adequate. So the profile can be 
a piecewise combination oi a^, crc, and crj^, except for 
scenario II, where only an applies. 
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In all cases, for the integral in Eq. (PT|) for (2(p) — 1) 
it turns out that at large L the contribution from gr 
dominates, and it alone gives the L — >■ oo value. This 
is because of the limitation of the range of the tan func- 
tion in ctl, to prevent its divergence, and of the range 
c^/^ cx L^/"^ ] of (T(7. This makes their contributions 
to the integral less than that from cr^ by a factor which 
vanishes as L increases. 

Note that, quite generally, the limitation of ctl requires 
xq to satisfy Xq — (— L/2) < 7r/q(— L/2); in the limit 
L ^ oo this restricts the variable A defined above to 
the two possibilities A = 1 (envelope apex very near the 
left boundary) or A > 1 (apex [well] outside). These 
are respectively Scenarios I and II, whose details are now 
exhibited. 



1. Scenario I 

In this case, where A = 1, we investigate its quantita- 
tive character and which values of (p) and J are con- 
sistent with it. Firstly, from Eq. A = 1 makes 
J — (1/8)(1 — 9). For possible values of (p), we consider 
Eqs. (PT|) and ([^^ . We chose the integration constant a 
so that a: = a is the center of the kink. If the kink is in- 
side the system, the further it is to the right the smaller 
will be the integral, and the associated {p) — 1/2. There 
is clearly a least (p) in Scenario I, applying when the kink 
is as far to the right as it can be (consistent with PBC). 
But Scenario II allows displacement of the envelope to 
the left (A > 1), and with fixed kink position this affects 
the value of the integral, since the more the envelope is 
displaced to the left, the larger will be the amplitude k{x) 
of the tanh at any particular x inside the system. 

So, small \{p) — 1/2| can be achieved with envelope 
apex near the left hand boundary, by adjusting the kink 
position (Scenario I), while (p) nearer or 1 needs a 
large displacement {0{L)) of the envelope to the left, 
corresponding to A > 1 (Scenario II). 

For illustration, consider the special case (p) = 1/2 for 
which the numerical profile is actually shown in Fig|T] 
In the Figure it is evident that the required zero value 
of the integral between the profile curve and the x-axis 
is achieved by having the abrupt rise of the curve, cor- 
responding to the kink, where it is. To the right(left) of 
the kink the curve follows the upper (lower) branch of 
the envelope function {k{x) is monotonically increasing). 
At the extreme left is the region around xq (necessarily 
small) where the tanh has become tan; its near divergence 
makes it easily able to match the PBC requirement. Thus 
one sees, in retrospect, that the fit shown in Fig (T] is in 
fact quite misleading. 

Scenario I is consistent as long as the kink stays within 
the right boundary of the system. Then, (t(L/2) at that 
boundary is positive and the PBC requiring ct{~L/2) to 
have the same positive value can be readily satisfied, as 
the tan form needs only a very small adjustment of its 
argument (within « tt) to achieve this. At the same time 
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Figure 3. (Color online) Kink position a/L against density 
(p) [see Eqs. (gZl) and pg)) ], for systems with PBC and rate- 
gradient values 9 as shown. 



the spatial range in which the tan form applies has to 
be very small to prevent unphysical a^s. Of course 
and (7(7, and the relationship of their constants to a are 
needed to complete the determination of the profile. 

This discussion is easily generalized and made more 
quantitative by using the integration result in Eq. (j22ll . 
with the appropriate real version K{x) of 7(x), together 
with the fact that when the kink aiX = a— xq = A 
lies well inside the system K{X) — K{A) is large (©(i)) 
at both limits, but of opposite signs. Further, the 
kink width (= w), such that the argument of tanh in 
Eq. (1^^ changes by 0(1) between x = a ± w/2, is 
w sa [dk{A)/dA]-^ = [(A -I- c) /A] 1/2, which is 0{l) for 
a = 0{L), except near xq where w diverges. Hence the 
integration result is, in the limit of large L, 

2((P)-^) = | (KiL)-2KiA)) , (27) 

where 

K{X) = [X{X + c)]i/2 _ ctanh^i 



X 



X + c 



1/2 



(28) 



For the special case (p) = 1/2 the kink position then 
has to be such that K{A) — (\/2)K{L) which, using 
the explicit form of K [see Eqs. (l26t and ((28)) ]. gives 
A/L = (1/2) + a/L « 0.61 for 9 = 0.2, consistent with 
the kink position in Figure [1] The general solution for 
the kink position against particle density is exhibited in 
Figure [3l for illustrative values of 9. Note that the range 
of values of (p) for which solutions are found is symmetric 
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relative to {p) = 1/2, and gets broader with increasing 
[see also Eq. below]. 

Larger values of a are associated, through Eq. (E^ . 
with (p) < 1/2, up to the limit a — L/2 when the kink 
center is at the right boundary. Then the mean density 
takes the limiting value {p)c such that 



{{P)c 



1 

2L 



K{L) 



r{[2(i 



l)]V2_(l_i)tanh-i 



29 



1 



1/2 



}. (29) 



This marks the condition where the two scenarios meet, 
and will correspond to the limit of a plateau region (in 
which, for L — oo, J — ^{1 ~ 9) applies) in the "funda- 
mental" diagram relating J with {p) and 9. 

At {p) > 1/2 one finds equations identical in form to 
Eq. (Ell), (EHl), and with \ - {p) replacing {p) - \. 
For small 9, Eq. ^ gives {p)^ - 1/2 = ^{^/2/i) 9^/^ . 
Thus, the extent of the plateau in the J — {p) diagram 
vanishes as — ^ 0. In this limit, for {p) still within the 
plateau, one can show that the height of the kink vanishes 
as 01/2. 



provides the following current-density relation outside of 
the plateau region 



1 



J = - (1 - A0) , 



(32) 



with 



iItt - <P» = [(-^ + iXt; + 1)1' - [(* - iKt; - 1)1' 



4 



A) {tanh 



" A + 1 ■ 


2 


■ A- 1 " 


— tanh~i 


e-i + 1 


0-1 - 1 



}• 

(33) 



A similar procedure applies for the complementary sub- 
case, (1 — {p)) < {p)c < 1/2, by particle-hole duality. 

Eqs. and ([55)1 can be combined to give J as a 
function of (p) in Scenario II, for fixed 9. The range 
of values of (p) for which physically acceptable solutions 
are found is complementary to that limited by Eq. (P^. 
which marks the extremes of validity of Scenario I. 



3. Weak-bond interpretation of plateau current 



2. Scenario II 

Scenario II applies at (p)'s so small (for a given 9) 
that the kink center is beyond the right boundary of the 
system (see, e.g., the curves for (p) — 0.25, 0.125 in Fig- 
ure H]) . Then the apex position xq of the tanh envelopes 
has to go outside of the system on the left, and there 
has to be a small upturn in a at the extreme right of the 
system to satisfy the PBC, so the start of the kink is just 
visible there in Figure [5J and the kink center is actually 
beyond the right boundary. This means that the tanh 
profile applies throughout the system: 

cr{x) =(JB. = k{x) tanh(X(.T) - K{a)) , (30) 

where K{x) is again as in Eqs. (IMl) and (|26p . We will 
now have xq < —L/2 and L/2 < a < L/2 + w, where w 
is the kink width (of order 1). 

As discussed above, for given 9 specifying (p) < (p)c < 
1/2 [see Eq. ] will lead to xq/L < -1/2, so making 
J < (1/8)(1 - 9) and A > 1. As before, we use Eqs. ^ 
and (P^ . But now, since xq < —L/2, and for all x in 
the system x < a, K{x) — K{a) is at both limits neg- 
ative (and large) . So we have [ ignoring contributions 
to 2((p) — 1/2) L of lower order in L (from corrections 
to the adiabatic approximation, and from width of the 
kink), and the comparable small distance the center lies 
beyond the right boundary]: 

2((P) ~\)^\{ + 1)) - K{^{\ - 1)) } (31) 

where K{X) is as in Eq. ([28)) . This gives A in terms of 
(p) and 9, for (p) less than the critical value, and hence 



Before moving to numerical results, we introduce an 
additional piece of mean-field theory which will be useful 
later. 

As remarked above, a plateau current, similar to that 
predicted in Scenario I, is found in the TASEP with ran- 
dom rates p{x) [1, [l^, |Tl|. There, an interpretation is 
given in terms of the current limitation provided by the 
weakest bonds, p^, which suggests that the "maximal" 
current satisfies Jmax < Pw/^- The following generaliza- 
tion provides a direct interpretation and confirmation of 
the result Jmax — — 6*) predicted for the plateau phase 
in Scenario I. 

In the continuum mean field formulation, Eqs. ([3]), (0]), 
and ([5]) give for all x: 

J=ip(x){l-a2(x)-|^}. (34) 

The most limiting rate, occurring at a; = —L/2, is — 
(1 — 0), so Jmax is obtained from applying Eq. ([M)) 
there. In Scenario I, with A = 1, the tan solution Eq. 
applies in that region [i.e. X — 0{\)], which yields for 
the right-hand side of Eq. (jM]), using Eqs. (UHl), dM]), and 
(ESI): 

\p^{l + q') = \p^.{l- = {1 + 0{j)} , 

(35) 

hence confirming the infinite-system maximal current 

D. Numerical results 

We considered lattices with L = 2™ sites, 8 < m < 13. 
A time step is defined as a set of L sequential update 
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attempts, each of these according to the fohowing rules: 
(1) select a site at random; (2) if the chosen site, here 
denoted by x, is occupied and its neighbor to the right is 
empty, then (3) move the particle with probability p{x). 
Thus, in the course of one time step, some sites may be 
selected more than once for examination, and some may 
not be examined at all. 

We have found that the time needed to attain steady- 
state flow varies roughly with L^^^, similarly to the 
uniform-rate case [llij, for which this is well-known { 28l . 
SjUHIi and is in agreement with the correspondence be- 
tween the (uniforni) TASEP and evolution of a KPZ in- 
terface [3l[3l|,|3i,[33|. 

For density profiles, local densities were usually aver- 
aged over snapshots (taken at appropriately long times) 
of 10^ independent samples. For example, for L — 256 we 
found that steady state has been reached by time i = 10'' 
in most cases, except for points on the coexistence line 
for open BCs (see Section ITVT) where the approach to sta- 
tionarity is markedly slower. Although finite-size effects 
can be observed, they are generally small and act towards 
making any kinks sharper, relative to system size, with- 
out any qualitative change. Thus we can be confident 
that no significant physical features are missed by gen- 
erally exhibiting profiles corresponding only to L = 256, 
as done here. 

Figure[4]shows steady-state density profiles for 9 = 0.2, 
which although still in the scaling regime is a relatively 
steep gradient. The behavior is in full agreement with 
the theory developed in Subsection III Bl (i) according to 
Scenario I, there is a common envelope, pinned to the left- 
hand extreme of the system, for intermediate densities 
roughly between 0.3 and 0.7; (ii) within this range of 
densities a kink is present, whose location varies against 
(p) as predicted by Eqs. ((27)) and ([28)) : (iii) for densities 
further removed from 1/2, Scenario II takes over, and 
profiles follow either the lower branch of the envelope 
(with its (p)- dependent displacement) with an incipient 
kink at the right boundary (for (p) < 1/2), or the upper 
branch, in this case with a narrow downward turn at the 
left edge in order to satisfy PBC {{p) > 1/2). 

Envelope functions are already familiar in the profiles 
of constant-rate asymmetric exclusion processes (e.g., on 
the coexistence line), but they only involve new scenar- 
ios when their delimitation of density profiles is space- 
dependent (as above or, e.g., in asymmetric exclusion 
problems with Langmuir dynamics [38]). 

There are slight numerical discrepancies between pre- 
dictions of Subsection IIIBI and the data displayed in 
Fig. 21 which exemplify the finite-size effects referred 
to above. For instance, according to Eqs. (|77)) . ([^5| . 
and (123) [see also Fig. [3|], for 9 = 0.2 Scenario I should 
hold for 0.293 ■ ■ ■ < {p) < 0.707 .... However, the pro- 
file for (p) ~ 5/16 — 0.3125 already shows some devia- 
tion from the common envelope. Overall, we have found 
that the quantification of finite-system corrections, to- 
gether with accurate analysis and extrapolation to the 
L ^ oo limit, can best be accomplished when dealing 
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Figure 4. (Color online) Steady-state local density profiles for 
system with L — 256 and PBC, 9 = 0.2. Curve labels denote 
average particle densities. 



with steady-state currents, as shown in the following. 

Evaluation of steady-state currents involved averaging 
over Ns — 100 independent samples, for each of which 
Nc — ttL L"^^^ successive instantaneous current values 
were accumulated. We took ~ 130 for L = 256 and 
512, and 13 for larger L. The instantaneous current is 
n-movcs/ L, where rimoves is the number of particles which 
undergo successful move attempts in the course of a unit 
time interval, i.e., L stochastic site probings as defined 
above. As is well known [39], the width 5, J of the dis- 
tribution thus found is essentially independent of as 

— 1/2 

long as Ns is not too small, and varies as N~^'\ With 
the parameters as specified here, we managed to keep 5 J 
well below the finite-size difference between J estimates 
for consecutive values of L (for fixed 9, p). The relevance 
of finite-size effects for currents is illustrated for p — 1/2 
in Fig. O where 9 is restricted to small values for clarity 
of presentation; one can see that the curvature present in 
finite-L data is essentially absent upon extrapolation to 
L oo. We now discuss guidelines for extrapolation of 
finite-system currents Jl to their thermodynamic-limit 
value Joo. 

In line with general finite-size scaling ideas, we at- 
tempted single-power fits of our sequences of finite-L cur- 
rent data with an adjustable finite-size scaling exponent 
-0, for all available pairs 9 and p. We denote by 9c{p) 
the gradient intensity value above which J{p, 9) becomes 
independent of p. 

So, 9 < 9c{p) corresponds to Scenario II of Subsec- 
tion [lirc] above, while 9 > 9c{p) is associated with Sce- 
nario I. Although 9c still carries an L-dependence (thus. 
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Figure 5. (Color online) Steady-state currents J against gra- 
dient intensity 6 for p = 1/2 and PBC, for system sizes as 
shown, plus extrapolated curve; for details of extrapolation, 
see text. 

e.g., the mergings of J — 6 curves shown m Fig. [2] take 
place at slightly different locations for L ^ 2048), it is 
a rather small effect compared to the overall range of 
0- variation investigated. 
Results were as follows: 

(1) For < 6* < Oaif), ~ 1 (Scenario II); 

(2) For > 0c{p), V ~ 1/2 (Scenario I). 

In the immediate vicinity of Oc{p), on both sides, we had 
rather serious convergence issues, so there we generally 
resorted to fixing tp — 1/2, for which the corresponding 
extrapolations fell in smoothly with the remaining ones 
outside that interval. For region (2), we estimate the 
uncertainty for -0 to be of order 10% at most. 

Thus, for the extrapolated points in Fig. [5j V' G 
(0.45, 0.55) was found in all cases except that correspond- 
ing to 6 — Q, for which -ij) k, 1. The case of {p) =1/2 
shown in that Figure is somewhat exceptional in that, as 
remarked at the end of Subsection IIII C II above, there 
the extent of validity of Scenario II corresponds only to 
the hmit 6* ^ 0. 

In Fig. El below we present the set of extrapolated cur- 
rents for p — 1/2, corresponding to < 9 < 0.2, to- 
gether with the mean-field prediction of a straight line 
Jmf{9) — I (1 — 6*) for Scenario I (see also the weak-bond 
interpretation given in Subsec. IlII C 3|) . The agreement 
is remarkable. 

Considering now the extrapolated currents for (p) ^ 
1/2, one sees in Fig. [7] that the variation of J against Q 
is generally much slower where Scenario II holds. In the 
vicinity of Bci^p)^ due to the convergence issues mentioned 
above, we considered systems of sizes up to L = 8192 




Figure 6. (Color online) Points are extrapolated currents for 
system with p — 1/2, PBC. Full line is the mean-field approx- 
imation (see text). 

(away from that region, we found that using L < 2048 
was generally enough to distinguish a reliably smooth 
trend as L^^ — > 0). Upon extrapolation we found the 
small overshoots shown in the Figure, which when trans- 
lated to J — (p) diagrams for fixed 9, would amount 
to reentrant behavior. For the largest deviation found, 
corresponding to p = 5/16 at 6 = 0.175, one gets 
J — 0.1039(1), while the value for p = 1/2 at the same 6 
is 0.1030(1). Although the average values differ by just 
under 1%, when converted in terms of (estimated) un- 
certainties this difference is equivalent to nine error bars. 
So, this effect appears to be real. 

The data in Fig. [7| can be used to test Eq. ((29)) . In 
order to do so, for fixed (p) < 1/2 one needs to establish 
the boundary between the ranges of validity of Scenarios 
I and II, as given by numerical simulations. Due to the 
overshoots just referred to, this task carries some am- 
biguity. For simplicity, we assumed such location to be 
where the respective J — 9 curve first crosses that for 
{p) — 1/2, upon increasing 9. Fitting the data thus ob- 
tained to the form ^ — {p) — a9^ , one finds a — 0.471(4), 
b = 0.51(1). These are to be compared, respectively, 
to a = V2/3 = 0.4714 ...,6=1/2, from the small-61 ex- 
pression of Eq. (|29p (see paragraph below that Equation) . 
Thus, the above assumption seems justified. 

Furthermore, data within the region of validity of Sce- 
nario II can be compared with the predictions of Eqs. (|32p 
and dMl)- We used 9 = 0.06. One sees in Fig. [5] that the 
agreement between theory and extrapolated numerical 
results is indeed excellent. The prediction of a plateau 
for Scenario I is also borne out by numerics, within er- 
ror bars. One cannot see unequivocal evidence here for 
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Figure 7. (Color online) Extrapolated steady-state currents 
versus gradient intensity for systems with PBC, and densities 
as shown. Note overshoots. 



a reentrant behavior similar to that found in Fig. [71 It 
is possible that such an effect, if present, is smaller than 
for the cases depicted in the latter Figure. This would 
be in line with the observation that the amplitude of the 
reentrance decreases with decreasing 9 . 

Fig. [5] also shows the J — {p) relation for uniform 
hopping-rate systems, for comparison. 

A current-density diagram very similar to Fig. [8] was 
obtained in Ref. Il4l . for the partially asymmetric exclu- 
sion problem with spatially-varying hopping rates. 



IV. OPEN BOUNDARY CONDITIONS 

A. Introduction 

With open boundary conditions, the following addi- 
tional quantities are introduced: the injection (attempt) 
rate a at the left end, and the ejection rate /3 at the right 
one. Calling p^, pn the stationary densities respectively 
at the left and right ends of the chain, one has for the 
current J at the boundaries, and anywhere inside: 



a{l - Pl) ^ J ^ P PR 



(36) 



Scenarios I and II, regarding the existence and location 
of an "envelope", discussed in the preceding section, still 
apply here, with similar consequences upon the system- 
wide current. The overall picture turns out to be rather 
like that for open systems with uniform hopping rate [if- 
0) l40l - l43l |: a maximal-current phase arises for suitably 
large a, /3 (where Scenario I takes hold); elsewhere, one 



Figure 8. (Color online) Points are extrapolated (L — > oo) 
steady-state currents; full line is J — (p) relationship from 
Eqs. (|32|) and ((33} [Scenario II]; long-dashed horizontal line 
is J = 0.1175, from Eq.((25ll with A = 1 [Scenario I[. Short- 
dashed line is J — {p) relationship for uniform hopping-rate 
systems, Eq.(l2|). 



has less-than-maximal current, although with either low 
or high density, the latter two subphases being separated 
by a coexistence line; Scenario II applies. See Fig. [HI and 
corresponding insets. 

The robustness of the three-phase structure in the 
present case is in line with the results of Ref. [ij. In 
their study of the partially asymmetric exclusion prob- 
lem, with spatially-varying right- and left- hopping rates 
p{x) and qix) respectively, those authors always found 
three phases, as long as p{x) — q{x) did not change sign. 

In the maximal current phase, the following specific 
features are noteworthy: 

(i) the steady-state system-wide density (p) is very close 
to the value which, for PBC, corresponds to the lower 
limit of validity of Scenario I. This is because, from the 
conditions giveir in Eq. (j36p . for large a,/3 one must 
have PL "large" and pn "small". Thus the density pro- 
file essentially follows the lower branch of the envelope 
function. Slight departures from that occur within short 
("healing") distances from the extremes, in order to com- 
ply with the exact values dictated by Eq. ([55)) . The latter 
effects account for the fact that {p) is not strictly con- 
stant throughout the maximal-current phase. Although 
Eq. ([36[) imposes the same constraints for systems with 
uniform hopping rates, there the envelo pe is trivially 
x—independent, and (p) is close to 1/2 l27l l4l|. 

(ii) In contrast to Scenario I with PBC, the steady-state 
profiles here do not show a kink inside the system. 

(iii) Similarly to Scenario I with PBC, the tan- like seg- 
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ment of the profile at the extreme left of the system is 
essential in the local density adjustment near that edge. 
However, as just mentioned, such adjustment is here im- 
posed by Eq. as opposed to the former case where 
the constraint arises from demanding continuity of p to 
obey PBC (combined with the existence of a kink further 
to the right). 



B. Theory and Scenarios 

In the "low current" Scenario II, with the apex left of 
the system's left boundary (xq — — AL/2, with A > 1), 
one has only tanh type solutions for all x\ a — 2{{p) — i) 
is limited by the envelope ±k{x): 



k{x) 



4J 



p{x) 



1/2 



{20/L){x--xo) 



1 + {26/L)x 



1/2 



(37) 



So, 



J=^(l-fc^(x)) = i(l + ^x) {l-k\x)). (38) 



Taking x = x^, where k(xa) = 0, gives 



J=-(l-A(?) , 



while evaluating Eq. ((57)) at x = gives 



(39) 



fc(T 



1/2 



(40) 



which sets the upper [p^ ^) and lower (p^ ^) bounds for 
the density at the extremes: 



(41) 



With a{x) — k{x) ta.nh{K{x) — K{a)), where a; = a is 
the position of the kink, and K — J k{x) dx, one can see 
from the insets in Fig. IHlthat, considering the situations 
corresponding to profiles types (i), (ii), and (iii) shown 
there, the following constraints hold: 



(i) 
(ii),(iii) 




(i),(ii) 
(iii) 



(42) 



So, using Eqs. p6p and (155)) we obtain, for the possible 
values of a and /3 in the three situations: 



J 



J 



Pi 



Pi 



-(l-0)(l-fci) 



1 



Pt 



%),(ii) = ^ = 7(l + ^)(l-M 
Pr ^ 



^(iii) e 



J J 



Pi 



[1 - kR, 1 + kl 



(43) 



For the coexistence line (CL), in which the kink lies 
wholly inside the system, i.e., profile type (ii) above, the 
results established in Eqs. (|43)) . together with Eqs. (|39)) 
and (HO)), give the current: 



JcL = a 1 - 



2a 



= /? 1- 



2/3 



1-ej ' \ 1 

and the equation for the CL shape as follows: 

2/32 



1 



1 



/3 



(44) 



(45) 



In this "low current" Scenario II, A > 1 and, from 
Eq. dMI), A < l/e. So the extent of the CL in {a.p) 
parameter space, and the current there, are limited to: 



< a 1 - 



2a 



= /3 1- 



2/3 

i + e 



= JcL < ^(1-^) 



(46) 

The same form of current, Eq. ([HH)) . and the same limi- 
tation 1 < A < 1/9, apply for the high density and low 
density sub-phases (corresponding to profiles of types (i), 
(iii)) which the coexistence line separates in this low- 
current Scenario II. Actually A = 1 is the boundary be- 
tween maximal (plateau) current phase (corresponding 
to profiles of type (iv)) and the lower current phase(s). 

Since for A = 1, fc^ = 0, fc^ = {29 / {I + e))^'"^ , using 
Eq. (|^5)) the phase boundaries are (in addition to the 
coexistence line): 

a = a,{9) , p > Pc{9) (47) 

(between subphase (iii) and maximal current phase), and 

13 = /3,(0) , a > ac{9) (48) 

(between subphase (i) and maximal current phase), 
where 



«c = i(l- 



29 



1 



(49) 



From Eqs. (|i5)) and (|1^ . the slope of the CL is unity at 
the origin, i.e. the same there as that for the uniform-rate 
case, and diverges at the endpoint (ac, /3c)- 
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(iv) 




Figure 9. Schematic phase diagram for TASEP with hopping- 
rate gradient for open boundary conditions. Locations of 
phase boundaries are 6- dependent [see Eqs. (|45|) - (|49[1 ]. CL 
stands for coexistence line (between high- and low- density 
phases). The insets show typical density profiles for each 
phase (see text). 

Given a and /3, Eqs. (gOl), and (gS]) give J and A, 
and then Eqs. (PT|) and can be used for the deter- 
mination of (p), anywhere on the phase diagram where 
Scenario II applies. For points on the CL, however, an 
adaptation is needed in order to account for the presence 
of a kink inside the system. Then, the amended form of 
Eq. ^ reads: 

2((P>-^) « \{K{^{\ + l)) + K{^{\-l))-2K{A)} , 

(50) 

where X — A \s the position of the kink. This can be 
found by keeping track of the leading finite-size correc- 
tions (from the asymptotic values ±1) to the tanh forms 
at the ends [ to obey the constraints given by Eq p6p ] . 
One gets the prediction {p) = ^ everywhere on the CL, 
for any 6. 

C. Numerical results 

In numerical work with open boundary conditions, we 
kept to 6* = 0.2. 

Initially we investigated the shape of steady-state pro- 
files deep inside the high-density, low-density, and maxi- 
mal current regions given in Fig. [51 as well as at a point 
on the CL at (a, /3) = (0.1, 0.087868) [about halfway be- 
tween the origin and the endpoint of the CL, see Eqs. (gS]) 
and (|49l) ] . We found profiles which conform respectively 
to types (i), (iii), (iv), and (ii) shown in the Figure, in 
agreement with the theoretical results given above. 

Deep inside the maximal-current phase, ai a — (3 — 
0.375, we carried out a finite-size scaling analysis of 
steady state currents, using systems with L < 8192. 



The finite-L values thus obtained were very close 
to those corresponding to PBC and = 0.2, for {p) in 
Scenario I. They approach the same extrapolated value 
Joo — |(1 — ^) found there, with the same type of finite- 
size corrections, i.e. Jl — Joo ~ L~'^ , -ip « 0.5. 

Elsewhere on the phase diagram, we calculated steady- 
state currents and densities at selected points, using only 
L = 1024. Results are shown in Table |T1 

For the point on the CL, the central value of (p) is 
close to 1/2, as predicted in Subsection lIVBi but the 
density fluctuations, associated to phase coexistence, are 
apparently very large. Related profiles (at the relatively 
small value 6* = 0.2 being used) are consistent with the 
kink being located with roughly equal probability any- 
where in the system, as in the 6 = Q case iEI. An 
approximate calculation, assuming this Q, and replacing 
the envelope by one with and Ur both set equal to 
their root-mean-square value, provides an estimate for 
the root-mean-square density deviation, (<5p)rnis = 0.17, 
in line with the result quoted in TablelD The relationship 
between the current JcL, Oi, and /3 given in Eq. (144p is 
verified to very good accuracy. 

In the high-(HD) and low-density (LD) phases, agree- 
ment between theory and numerics is excellent, in part 
because finite-size effects are small there, where Scenario 
II holds. 

For the set of three points inside the maximal cur- 
rent (MC) phase, the currents are indeed close to each 
other, their value differing from the infinite-system one 
^oo = i(l ~ ^) t>y well-understood finite-size corrections. 
The corresponding densities are also very close, and in 
good accord with the prediction that the corresponding 
profiles should essentially coincide with the lower branch 
of the envelope function. Recall that, for PBC, this is ex- 
pected to happen, at 6* = 0.2, for (p) close to 0.3 (see Sec- 
tion |TTT|. The larger spread between densities in the MC 
phase, when compared to that between currents, is to be 
expected (see comments in Subsection IIV A\ . One gets a 
smaller difference between numerical results and theoret- 
ical predictions by looking at a point on the borderline 
between LD and MC phases (LD/MC). Even then, the 
agreement is not as close as that found deep inside the 
HD and LD phases. Such effects reflect the cor- 
rections pertaining to Scenario I. 



V. DISCUSSION AND CONCLUSIONS 

We have developed a mean-field/adiabatic theory for 
the one-dimensional TASEP with smoothly- varying hop- 
ping rates. Its application to the uniform-gradient case 
is shown, upon comparison with extrapolations to the 



It is the number of particles present in the system, rather than 
the kink location, which is expected to have uniform distribution, 
so our procedure is not expected to apply at large 9. 
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Table I. Average steady-state currents J and densities (p) 
for systems with 6 — 0.2, L — 1024, and assorted injec- 
tion/ejection rates {a,/3). Jth and (p)th refer, respectively, 
to currents and average densities calculated by the theory 
given in Subsec. II V Bl Phases specified in column 1 are, re- 
spectively: CL: coexistence line; HD: high-density; LD: low- 
density; MC: maximal-current. Refer to Fig. |9] and text. 



Type a P J 

CL 0.100 0.087868 0.0750(3) 
HD 0.400 0.100 0.0834(4) 



LD 0.100 
LD/MC 0.200 
MC 0.450 
MC 0.650 
MC 0.250 



0.400 
0.400 
0.450 
0.250 
0.650 



0.0751(3) 
0.1015(2) 
0.1036(1) 
0.1036(1) 
0.1031(1) 



Jth 
_3_ 
40 
J_ 
12 
_3_ 
40 
J_ 
10 
J_ 
10 
J_ 
10 

J_ 

10 



(p> (p>th 

0.51(10) i 

0.781(1) 0.78241. 

0.188(1) 0.18839. 

0.297(1) 0.29245. 

0.3113(7) 0.29245. 

0.3120(6) 0.29245. 

0.3074(7) 0.29245. 



L — > oo limit of numerical simulation data, to give very 
accurate results. Evidence for this is exhibited espe- 
cially in Figs. H H [71 [i and Table [H Thus, for PBC 
it appears that the J — (p) — 9 relationship given by 
Eqs. ([5^ and ([55)1 is exact for Scenario II of a (p)- de- 
pendent current. While simulations essentially find the 
constant-current plateau predicted for Scenario I with 
PBC (at values of J in full accord with theory), a small 
amount of nonmonotonic dependence of J on (p), near 
the edge of the corresponding region, appears to be 
present. Although extrapolation of finite-system current 
results turns out to be plagued with convergence issues 
precisely in this region, a systematic trend is found to- 
wards increasing values of the calculated overshoot as 9 
increases (see Fig. [7]). Thus one cannot definitely discard 
the possibility that such overhangs are real effects. 

Being mean field in character, the theory presented 
here cannot predict, e.g., current fluctuations (i^ . [isj . 
nor fluctuation-related finite-size corrections. However, 
our numerical evidence shows that in the plateau region, 
i.e., within Scenario I (both for PBC and open BC's), the 
dominant finite-size current corrections are of order L^"^ , 
Tp « 0.5. This indicates that an additional mechanism is 
present, whose effects obscure the usual (uniform- hop- 
ping rate) fluctuation- induced L^^ terms [the latter are 
clearly identified in our numerics, not only for f? = 0, but 
also wherever Scenario II holds ] . 

The mean-field theory explains the L^^/^ corrections 
as arising from a tan —like part of the profile which 
lies inside the system only in Scenario I. This occurs 
near the envelope apex, in a region of width T:/q(X), 
where q = dQ/dX, see Eqs. dH]) and There, 
X = C(l), c = 0{L), hence from Eqs. 
and (HH), TT/q{X) = ©(L^/^). One must quantify the 
subdominant size-dependent effects originating from this 
region. 

For PBC, notice that in the range of x where the tan- 



like profile holds, a = 0(1) so it gives a contribution to 
the integral for L (p) of order L^/^, out of a total of order 
L. This provides a correction of relative size L^^/^ in (p) 
for given J. By inverting the J — (p) relationship (since 
for PBC it is the density which is fixed), one is left with 
the observed current corrections 0{L^^^^). 

The argument for open boundary conditions is slightly 
different, because (p) is not fixed by initial conditions 
and J is determined by the boundary injection/ejection 
rates. So, we look directly at the current and its rela- 
tionship with a and p^, as given in Eq. (j36p . since the 
tan solution applies near the left edge. By Eq. ([T^ . this 
is tr ~ — tan q^ix — const.), where q^ is the value of 
q = {-X/{X + c)y/^ ntX ^x~xo^ 0(1), with X <0. 
So qyj = 0(L-i/2), again by Eqs. and 
To provide the required injection current, one must have 
(T = 0(1) near the left edge, while only a small change 
Ax in position, of order Ax w 1/qw = 0{L^^'^) will take 
the tan — > 0. The upshot is that the average change in 
caused by a change of 0(1) in xq is ^ L^^/^. Hence with 
open boundary conditions, whenever Scenario I applies, 
the finite-size correction in J = ^ (1 — gl) is 0{L~^/^) . 

By similar arguments one finds that, when a kink is 
present, its width generally gives corrections of order 
to (p). On the other hand, corrections coming from the 
region where the validity of the adiabatic approximation 
breaks down are of order L"^/"^. Since these only oc- 
cur when the apex is inside the system, i.e., when the 
L^^/^ tan- originated terms are present as well, they are 
dominated by the latter. 

In closing, we note that a number of extensions of 
this study suggest themselves. Steady state behavior 
for other spatial dependences of rates, particularly wells, 
should be amenable to similar procedures. The same is 
true for studies of the dynamics. To develop the the- 
ory beyond the mean field limit is a more formidable 
challenge, but for slowly varying rates the adiabatic ap- 
proach should still apply, possibly combined with existing 
exact methods for uniform systems. Phenomenological 
domain- wall approaches (46l . l47l | would be a likely way 
forward. 
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